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Abstract. By employing the adaptive network simulation method, we demonstrate that the ensemble- 
averaged stress caused by a local force for packings of frictionless rigid beads is concentrated along rays 
whose slope is consistent with unity: forces propagate along lines at 45 degrees to the horizontal or vertical. 
This slope is shown to be independent of polydispersity or the degree to which the system is sheared. 
Further confirmation of this result comes from fitting the components of the stress tensor to the null stress 
'constitutive equation.' The magnitude of the response is also shown to fall off with the -1/2 power of 
distance. We argue that our findings are a natural consequence of a system that preserves its volume under 
small perturbations. 

PACS. 45.70.Cc Static sandpiles; granular compaction - 83.80.Fg Granular solids; granular compaction 

1 Introduction in a constitutive equation, as in elastic systems. Further- 
more, experiments measuring the stress under conical and 

It is becoming increasingly evident that the transmission wedge-shaped sandpiles have indicated a pressure mini- 

of stress through a disordered body of incompressible, co- mum below the a P ex when the material is poured from a 

hesionless beads obeys principles that are distinct from P oint ( or line ) source ' but not when [t is sprinkled uni- 

those of either elastic or elastoplastic materials Qggg§§ ■ formly over the base 1 ' This de P endence on construction 

Even defining a strain field for a granular body is prob- history is also quite unlike elastic or elastoplastic materi- 

lematic let alone relating the strain to the stress a ^ s - 
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To circumvent the difficulties arising from an ill-defined ulations for additional data. However, conventional simu- 

strain, a class of 'constitutive equations' have been pro- lation methods cannot directly probe the linear response 

posed that relate the independent components of the stress regime and must instead apply a finite force, which runs 

tensor (jy without reference to a displacement field (|^,^| . the risk of perturbing the system. To this end, two of 

(Note that, like many authors, we use the phrase 'constitu- us (AT, TW) have devised a novel simulation method, 

tive equation' even in the absence of a well-defined strain called the 'adaptive network algorithm,' which directly 

field or history independence). In their most general form, probes the linear regime in a system of frictionless and 

these equations suppose that the cry can be related in a rigid discs jllj . 

linear expression, with coefficients that in general depend T , , , , , n 

1 ' ° r In tms paper, we employ the adaptive network algo- 

on the construction history of the pile. The basic form ., , . . , , r , • c ■ 

J 1 ntnm to investigate tne stress response function m iric- 

of these equations has also been confirmed by mean field , ,. , , , , . , ,. n 

J tionless disc systems, and to ascertain it it obeys wave-like 

analysis of frictionless spheres, where it was referred to as , . TTr , . , , . a , ., . , , , , . 

J propagation. We begin by briefly describing tne simulation 

method in Sec. ^J, before demonstrating in Sec. [| that the 



the 'null stress' law 10 



A remarkable feature of the proposed equations is that null-stress law is obeyed, with the stress tensor approxi- 

they are hyperbolic in form, and thus belong to the same mately obeying the relation a yy — c 2 a xx with c « 1. This 

class of equations as the wave equation. Hence they pre- result has already been found [ pT| for weakly polydisperse 

diet that the response to a localised perturbation should systems, with the same value of c; here we extend the 

propagate outwards from the source of the disturbance result to strong polydispersity. 

in waves, with the direction of increasing depth playing mi <• , • •, u - i j • o HI i 

' -r i/ a 1 he response function itself is analysed m Sec. where 

the role of time. To test this prediction, experiments have ., . , , , , , ,,. , , , , , 

r it is shown to be concentrated along two light rays that 

recently been performed on systems of photoelastic discs , r e , , 

J propagate outwards from tne source of tne disturbance. 

that have a small load applied to the upper surface, and y-, ,, , j.,. , -, . r n . 

furthermore, the same speed of light cw lis found m all 

the stress response in the bulk measured HI . These seemed . ,. ,. r ,, £ 

r 1-J J cases, irrespective ol the direction ol the perturbing lorcc, 

to show wave-like propagation over all distances for low , ., e , , , , . , 

° tne polydispersity of the system, or tne degree to wnicn 

polydispersity, but only for short distances for high poly- , , . , . , , m, . , , . . , 

f j r ji j the system is being sheared. 1 his behaviour persists to the 

dispersity. To further complicate the issue, similar exper- , , , , , , , , . . „ 

^ J * it- largest depths we have been able to simulate, typically 

iments on sand have failed to show any trace of wave-like , lr . , r , , . .. , ,. 

around 1U — 15 bead diameters m tne vertical direction, 

propagation whatsoever Ml. u- u • L , .,, , , , ,. 

° 10 wnicn is comparable with the experiments on photoelastic 

Until more experiments have been performed and a discs; however, they found no evidence of c w 1. In Sec. |^ 

consensus reached, it is natural to turn to computer sim- we attempt to explain this discrepancy by arguing that 
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c — I corresponds to a system that conserves its volume 
under small perturbations. Finally, in Sec. |^ we summarise 
our results and suggest directions for future work. 

2 Overview of the simulation method 

The full derivation of the adaptive network algorithm and 
its associated assumptions has already been presented at 
length elsewhere pH ; here we only mention those details 
pertinent to the current work. The algorithm consists of 
two stages: firstly, a packing of discs is constructed, and 
secondly, the contact force network is generated in re- 
sponse to an applied load. We consider two-dimensional 
systems of frictionless discs a, with radii R a that are uni- 
formly distributed over the interval [1, i? max ]- The param- 
eter i?max > 1 will be used as a measure of the polydis- 
persity of the system. The pile is built upon a fixed rough 
base, consisting of a row of discs belonging to the same 
size distribution as those deposited, and whose centres all 
have the same vertical coordinate. Periodic boundary con- 
ditions are obeyed in the horizontal direction. For compu- 
tational efficiency, the beads are assumed to be weightless 
except for those at the free surface, which are subjected 
to a uniform load F cxt . Here, F cxt can point in a direction 
other than vertically downwards, thus allowing the system 
to be sheared. 

The pack of discs is built by adding each successive 
disc at the lowest available point. Once added, each disc's 
centre remains fixed throughout the subsequent relaxation 
process. The sequential nature of the initial packing means 
that it is straightforward to calculate the contact forces 
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as the applied load propagates from the upper surface to 
the base; see jy] for details on how this is done. In gen- 
eral, some of the contact forces will be tensile and thus 
incompatible with a cohesionless medium. Therefore we 
now modify the contact network until all tensile bonds 
have been removed. This relaxation process exploits the 
fact that only a fraction of the physical neighbours of a 
given disc are in contact with it. Stability is achieved by 
successively removing contacts with some neighbours and 
adding them to others. 

In a pack of discs such contact re-arrangements re- 
quire motions of the discs. To avoid such motion we make 
a restrictive geometric Ansatz. We suppose that the beads 
once assembled have been deformed so that all non-con- 
tacting neighbours are nearly in contact. Specifically, we 
add material to each pair of non-contacting neighbours so 
that the gap between them is infinitesimal. The remain- 
ing gap is a) at the centre of the original gap, b) oriented 
the same as the original gap and c) proportional to the 
width of the original gap. Under this Ansatz contact rear- 
rangements can occur with arbitrarily little motion. Thus 
the positions of the beads and the contact angles can be 
taken as fixed during the relaxation process. The Ansatz 
thus simplifies the relaxation behaviour greatly. It also 
abstracts the essential process of contact rearrangement 
from the inessential process of bead motion. 

The adaptive network approach has a number of ad- 
vantages over traditional simulation methods. It operates 
directly in the limit of infinitely rigid beads, and also in 
the isostatic limit . That is, the number of contacts is 
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just sufficient to determine the positions of the beads and weak compressive contacts that are proportional to the 

the magnitudes of the contact forces. In two dimensions weight of a single bead. Physically, these beads occupy 

an isostatic pack has on average 4 contacts per bead. A a position in the force network that is shielded from the 

further advantage of using this method is that the linear applied load by an 'arch' above the bead, and thus play no 

response regime can be probed directly, a feature that will role in the macroscopic propagation of the stress through 

be exploited in Sec. ^ of the current work. However, the the system. In the language of [ ]17[ | , they are spectator 

computational resources required do not scale favourably beads. The proportion of such beads is typically small but 

with the number of beads N, with the memory needed finite, roughly 2% for a low polydispersity i? max = 1.1, 

varying as 0(N 2 ) and the simulation time increasing as increasing to around 5% for i? max = 3. There is also a 

0(7V 4 ). Thus its benefits are restricted to somewhat small slight increase in their numbers for sheared systems, 
piles, typically N < 500 here. 

A sample output from our simulations is given in Fig. [l], 

3 The null-stress law 

which shows a highly polydisperse system with i? m ax = 3 
that is first stabilised under a vertical load, and then 

As mentioned in the introduction, there has been con- 

restabilised under a shearing force angled at 20° to the 

_— . siderable recent interest in deriving an equation relat- 

vertical |13f. The positions of the beads is preserved un- 

ing the components of the stress tensor for a granular 

der the change of load; this is the hallmark of the adap- 

body [p]j8[P,[l^,p^|Jl9|] . In this section we test a proposed 



tive network algorithm. However, the network of contact 

class of equations against our simulations. Similar results 

forces adjusts to support the new load, and moreover the 

have already been published in pd[ , but only for low poly- 
strongest contacts have a clear tendency to align with the 

dispersity. Here we extend the analysis to highly polydis- 

direction of F . No such alignment is apparent for the 

perse systems with i? max = 3, thus eliminating concerns 

weaker bonds. The appearance of concentrated chains of 

over the possible effects of crystallinity. 

force that align in such a way as to support the load has 

n , , , • • , , • , . . irniFTcirrHi The stress tensor an generally has 3 independent corn- 

also been observed m experiments and simulations p4 lolffl. JO J 



ponents in two dimensions (the absence of local torque 

Before presenting our quantitative results, we remark 

fixes o~ X y — lyx). Since the forces must balance at every 

that although the contact forces must be non-tensile, they 

point in the system, cry obeys the stress continuity equa- 

need not be compressive. Highlighted in Fig. [j] are a small 

tions 

number of contacts with exactly zero force. It should be 
recalled that the beads in the bulk are weightless. In a 

real packing, these zero-force bonds would be replaced by c^ov,- = F° xt . (1) 
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This gives two equations in three unknowns, so an ex- large angle load is still linear, as indeed it must always be 

tra equation is required for closure. One proposal for this for any packing constructed by the adaptive network algo- 

'missing' equation, called the 'oriented stress linearity' rithm) . Equation (|^) may hold even for large angle shears 

model [|l| or 'null stress' law [[To), is that there exists a if is some function of the ratio F® xt /F| xt , but our data 

linear relationship between the u^, is too noisy to suggest a meaningful fit. 

Additionally, the finding that 77 w 1 is in agreement 

a ' xx = r/ + /i ° xv . (2) with earlier work on weakly polydisperse systems Jj"l| ]. 

Oyy Oyy 

This suggests that the 'speed of light' for stress propa- 

The coefficients 77 and fi depend on the construction his- 

gation is approximately 1, irrespective of polydispersity. 

tory of the pile, so /1 will vanish if the construction respects 

The analysis of the response function presented in the 

left-right symmetry, for example. 

next section also agrees with this observation. We post- 
To test the proposed equation (Q), we plot in Fig. ^| 

pone further discussion of this point until Sec. ||. 

the components of <7ij for systems under differing degrees 
of shear. Here the cry have been averaged over the bulk 

of the system. This is a valid procedure because, with our ^ ^^ponse function 

geometry and boundary conditions, the try are constant The previous section consi dered only the coarse-grained 
everywhere except the narrow surface layer. It is immedi- componcnts of a .. averaged over the whole bulk of the 
ately apparent from the data that a xy ja m = /F% packing. To gain further insight into the nature of the 
to very good precision. This is a trivial result that is eas- stresg propa gation, we now turn to consider the response 
ily derived from (@), and can be regarded as a consistency of the gystem to an mnnit esimal force F 7 applied to a 
check for our simulations. single bead 7 

More revealing is the data for normal components of xhe pre diction of the equation (§) for this problem is 
the stress, which approximately obey a xx = a yy for un- already well documented (see e.g. §0), but in brief, the 
sheared or slightly sheared systems. This suggests that substitution of @ into the stress continuity equations (@) 
H = 0, confirming that our construction procedure pre- gives 
serves left-right symmetry. The linear relationship between 
a xx and a vy breaks down under higher angles of the ap- 

dxiwyy + V>Vxy) + d y a xy = F° xt , (3) 

plied load, indicating that the shear has changed the struc- 

. , d x a xy + dyOyy = F^ xt . (4) 

ture of the system m such a way that /i is no longer zero, 

and left-right symmetry is violated. (Note that the re- After some straightforward manipulations, these can be 
sponsc within a system that has been stabilised under a rewritten as 
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(d y - c + d x )(d y - c~d x )(Tij = , (5) 

where = |(— fx ± \//i 2 + 4ry). These equations are hy- 
perbolic, and as such belong to the same class of equa- 
tions as the wave equation. Thus they predict that the re- 
sponse to the perturbation should propagate downwards 
only along the two characteristics, or 'light rays,' defined 
by Ay = c^Ax, where (Ax, Ay) is the relative displace- 
ment from the point of the perturbing force. 

If this picture is correct, then only those beads a and 
j3 with a point of contact that lies near to one of the light 
rays can have their contact force altered by the application 
of the perturbation. To test this prediction, we define the 
response function G(a[3\ r y) by 

Af aP = G(o/3| 7 ) • F 7 , (6) 

where Af a p is the change in the magnitude of the con- 
tact force between beads a and f3 due to the perturbing 
force F 7 acting on bead 7 . Thus G(a/3| 7 ) encodes the re- 
sponse to the force F 7 without requiring the direction of 
F 7 to be specified. Note that for our geometry-preserving 
procedure, only the magnitude of the contact force can be 
varied, as the angles were fixed during the initial prepara- 
tion stage. 

G(a/3\j) is plotted in Fig. || for systems of three dif- 
ferent polydispersities that have been stabilised under a 
vertical load. For the lowest polydispersity i? max = 1.1 the 
response is concentrated along rays radiating at f=s 45° an- 
gles, suggesting that c w ±1. This cannot be due to crys- 
tallinity, as the contacts in a close packed two-dimensional 
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packing lie rather at a 30° angle. For higher polydisper- 
sities, the light rays become less distinct due to the in- 
creased noise, and the magnitude of G(aj3\'y) decays more 
rapidly away from bead 7. However, on length scales of ^5 
beads diameters or more it is still possible to discern max- 
ima in G(a/3| 7 ) along the same 45° lines as before. 

To make these observations more quantitative, we plot 
m Fig. I the magnitude of the response |G(Z\x, Ay) | at dif- 
ferent depths Ay below the perturbation. It is evident that 
the peak response lies near to the lines Ay = ±Ax, again 
confirming that c w ±1. Furthermore, the amplitude of 
the peaks decays like (Ay)^ 1 / 2 , which presumably arises 
from the diffusive spreading-out of |G(Z\x, Ay)\ with in- 
creasing depth ]2ti|] . Diffusive spreading should also mean 
that the width of the peaks increases like (Ay) 1 / 2 , but our 
data is too noisy to test this claim at present. 

It is interesting to compare our results to the recent 
model of Bouchaud et al. which allows the force to prop- 
agate along straight rays until they are scattered by a 
random distribution of defects [|lj . These authors found 
that disorder causes a crossover from a two-peak response 
function on short length scales to a single-peak, 'pseudo- 
elastic' behaviour on longer scales. By contrast, we have 
found no evidence of any crossover, and the two-peak form 
of |G(Ar, Ay)\ persists to the longest length scales we 
have been able to simulate. It is not clear if this discrep- 
ancy has anything to do with the differing underlying as- 
sumptions between the two models, or if we are simply 
unable to generate sufficiently large systems to observe 
the crossover. 
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The response function for systems that have been sta- 
bilised under a shearing load are given in Figs. [| and [j[ 
It is readily apparent that applying the shear does not 
alter the orientation of the light cone, i.e. the characteris- 
tics have the same gradients rj ±1 as in the unsheared 
system. However, the magnitude of the response becomes 
stronger along the ray that extends in the same direction 
as the shear, and weaker along the opposite-pointing ray, 
as required by (||) . With increased polydispersity, there is 
a greater degree of noise and the weaker ray is eventually 
lost in the background, as seen in e.g. Fig. ^|c. 

5 Discussion 

The central finding of this paper has been the repeated ob- 
servation of a 'speed of light' c = 1 for the characteristics 
along which the stress propagates. This was first inferred 
from the form of the 'constitutive equation' described in 
Sec. ^j, which for unsheared or mildly sheared systems were 
consistent with a xx — rja yy with i] = c 2 rj 1. Confirma- 
tion came from the response functions in Sec. ^, which 
were concentrated within a 'light cone' with a central axis 
running parallel to gravity. The gradient of the surfaces 
of this cone were again « ±1. In other words, the two 
characteristics of the hyperbolic equation for stress are 
nearly orthogonal. Our earlier study indicates that 
this property survives even when the preparation proce- 
dure is not mirror-symmetric. In that case, the null-stress 
law has the general form a xx — rjcr yy +fj,(7 xy , with the value 
of rj still consistent with 1 . As seem from (|J) , if rj = 1 then 
the characteristics obey c + c~ = — 1 for all fj,, i.e. they are 
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orthogonal, just as in the fixed principle axes theorem 
The fact that the orthogonality of the characteristics ap- 
pears to be so robust, demands an explanation. 

Equality between the normal stresses a xx — a yy is 
also found in liquids, where it reflects the isotropic na- 
ture of the medium. However, this interpretation cannot 
be applied to our model, since the construction history 
and boundary conditions give rise to a system that is any- 
thing but isotropic. The anisotropy is most explicit in the 
response functions plotted in e.g. Fig. ||, which clearly 
show that the orientation of the light cone is fixed with 
respect to gravity, irrespective of the direction of the per- 
turbing force F 7 or the degree to which the system is being 
sheared. Clearly we must look elsewhere for an explana- 
tion. 

In the absence of a more compelling reason why c = 1, 
we tentatively propose the following explanation. Firstly, 
observe that the relation a xx = a yy means that the vol- 
ume of the system is conserved under small perturbations. 
To see this, consider a system that is contained within a 
rectangular box of width X and height Y. Let us suppose 
that, under the action of an external force f x applied to 
one of the side walls, the width of the system changes to 
X + 5X. Similarly, a force f y is applied to the top of the 
box, and the height changes to Y + SY. If both 5X and 5Y 
are sufficiently small, the topology of the contact network 
is not altered during this motion and the perturbation is a 
'soft mode' p0| , meaning that it does not change the total 
energy of the system and hence has no restoring force. This 
follows from the isostatic nature of the network; see [[To) 
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for details. To avoid runaway motion, the total work done 
on the system must vanish, i.e. f x SX + f y 5Y = 0. But the 
relation a xx = a yy can be rewritten as f x /Y = f y /X, giv- 
ing YSX + XSY = 0. The left hand side of this equation is 
just the change in volume of the system, which therefore 
remains fixed, as claimed. 

It is plausible to suppose that our observation of c = 1 
arises from a state of minimal volume or maximal com- 
paction, as suggested above. This suggests that c may be 
unity in a broader class of physical systems. Still, we have 
not demonstrated the necessary conditions for a minimal- 
volume packing. It may require a random close-packed 
initial state. It also may require frictionless contacts or 
regular bead shapes. Further, we have not demonstrated 
that our simulated bead pack has minimal volume in the 
requisite sense. Nor have we shown that it would continue 
to have c = 1 if we allowed the beads to move during 
relaxation. Thus the class of systems with c = 1 beyond 
that presented here cannot be fully characterised at this 
time. 

6 Summary 

In summary, we have investigated the propagation of stress 
through a body of frictionless rigid discs using the adap- 
tive network algorithm. This algorithm first generates a 
two-dimensional packing via a sequential deposition pro- 
cedure, and then calculates the contact forces in response 
to a load applied to the free surface. Tensile bonds are 
removed by altering the topology of the contact network, 
in a manner that preserves the mean coordination num- 
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ber; the packing is always isostatic. This is performed in 
the linear regime in which bead deformation and motion 
is neglected, so that although the topology of the network 
can evolve, the geometry is preserved. 

The linear response to a force applied to a single bead 
was calculated, and was found to propagate in two downward- 
pointing 'rays,' Our central result is that the 'speed of 
light' c of these rays is approximately 1, irrespective of the 
polydispersity of the bead size distribution, or the degree 
to which the system is being sheared. We have also argued 
that c = 1 may be understandable if the system naturally 
adopts a local minimum of volume. Furthermore, we have 
shown that the magnitude of the response function decays 
in a way that suggests a diffusive spreading-out with the 
vertical distance from the source of the perturbation. It 
would be interesting to see if other microscopic models 
also scale in a similar manner. 

There are many possible directions in which the re- 
sults presented in this paper can be extended. In partic- 
ular, we feel that there would be considerable advantages 
to devising a modified adaptive network algorithm that 
scales more favourably with system size, whilst still di- 
rectly probing the linear response regime. This would help 
to reduce the significant error bars found in some of our 
data. One way to achieve this might be to change the rule 
that updates the system when a single bond is changed, 
which is currently global and requires 0{N 2 ) calculations. 
Replacing this with a local rule should significantly im- 
prove the simulation times and hence the statistics, al- 
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though it is not clear what approximations such an opti- 4. R. Brockbank, J. M. Huntley and R. C. Ball, J. Phys. II 7 

misation would entail. 1521 (1997). 

The current work has focused purely on the problem 5. L. Vanel, P. Claudin, J.-P. Bouchaud, M. E. Cates, E 

of stress propagation through granular materials. How- Clement and J. P. Wittmer, Phys. Rev. Lett. 84, 1439 (2000) 

ever, granular media make up only a small subclass of 6 ' J ' Gen §' D ' Howe11 ' E ' Lon S hi ' R ' R Beh ™ger, G. Rey- 

, , .11 <• i , i ,i dellet, L. Vanel, E. Clement and S. Luding, preprint ]pond- 

systems that can be described as jammed, a class that c 

also includes glasses and foams p^, for example. It has 



mat/001 21 gj, 



7. G. Reydellet and E. Clement, Phys. Rev. Lett. 86, 3308 

recently been demonstrated that the distribution of inter- 

(2001). 

bead forces in these systems exhibits a kind of universality 

8. J. P. Wittmer, M. E. Cates and P. Claudin, J. Phys. I 7, 39 

when in their jammed state [p2| . Clearly the adaptive net- 

(1997). 

work algorithm could be a useful probe of the extent and 

9. J. P. Wittmer, P. Claudin, M. E. Cates and J.-P. Bouchaud, 

limits of any such universality. Since this is a separate 

Nature 382, 336 (1996). 

question to that considered in this paper, we will present 

10. A. V. Tkachenko and T. A. Witten, Phys. Rev. E 60, 687 

our results on this issue elsewhere |E3f. 
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(a) 



V V V V V V V 




(b) 

Fig. 1. Sample output of the adaptive network algorithm for 
a system subjected to (a) a vertical load, and (b) for the same 
packing under a shear of approximately 20°. The arrows point 
in the direction of the loading force F oxt . The black lines be- 
tween beads denote compressive contacts, with a force whose 
magnitude is proportional to the thickness of the line. The 
thick gray lines represent tenuous contacts of exactly zero force. 
Note that contacts between slightly separated beads are al- 
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Fig. 2. The components of the stress tensor as a function 
of the applied shear for an N = 400 bead system of polydisper- 
sity .Rmax = 3. The error bars for the a xx /a yy points are the 
standard deviation over 25 different packings. For comparison, 
the dashed line represents a xx — a yy . The o xy /a yy data lie on 
a line of slope 1 and the errors are smaller than the symbols. 
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(c) 



Fig. 3. The coarse-grained response function G(a/3\j) for 
packings under a vertical load and with polydispersity 
(a) i? max = 1.1, (b) i? max = 1.5 and (c) i? max = 3. Here, 
G (a/3 1 7) has been averaged over all those beads 7 whose cen- 
tres lie in a narrow horizontal strip approximately two thirds 
of the way up the pile. This central bead is represented by the 
large grey filled circle. The vectors shown represent G(a/3\'y) 
as a function of the relative displacement (Ax, Ay) from the 
centre of bead 7 to the point of contact between beads a <-> f3. 
Each vector points away from the small grey circle. The data 
was averaged over 100 runs of N = 500 beads. Each mesh point 
is separated by 1.5 bead diameters. 
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Fig. 5. The response function G(a(3\~/) when the applied load 
is angled at 20° to the right of the vertical, for polydispersities 
(a) -R m ax = 1.1, (b) iCax = 1-5 and (c) i? max = 3. 



Fig. 4. The magnitude of the response function |G(a,/?|7)| 
scaled up by a factor (Z\j/) 1//2 , plotted against Ax/ Ay, where 
{Ax, Ay) is the relative displacement from the source of the 
perturbation measured in units of the mean bead radius. Pos- 
itive Ay correspond to points below the source. Each line rep- 
resents the average over a horizontal strip of width 3, with a 
mean Ay as indicated in the legend. The polydispersity was 
(a) .Rmax = 1.1 and (b) -R max = 1.5, and N = 500 in both 
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Fig. 6. The same as Fig. ^| for a system under a shearing 
load, angled at 20° to the vertical. Only the results for weakly 
polydisperse system i? max = 1-1 is given, as the error bars were 
much larger than the data points for higher polydispersities. 



